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We analyze a system of two colliding ultracold atoms under strong harmonic confinement from 
the viewpoint of quantum defect theory and formulate a generalized self-consistent method for de- 
termining the allowed energies. We also present two highly efficient computational methods for 
determining the bound state energies and eigenfunctions of such systems. The perturbed harmonic 
oscillator problem is characterized by a long asymptotic region beyond the effective range of the 
interatomic potential. The first method, which is based on quantum defect theory and is an adap- 
tation of a technique developed by one of the authors (GP) for highly excited states in a modified 
Coulomb potential, is very efficient for integrating through this outer region. The second method is 
a direct numerical solution of the radial Schrodinger equation using a discrete variable representa- 
tion of the kinetic energy operator and a scaled radial coordinate grid. The methods are applied to 
the case of trapped spin-polarized metastable helium atoms. The calculated eigenvalues agree very 
closely for the two methods, and with the eigenvalues computed using the generalized self-consistent 
method. 



I. INTRODUCTION 



An understanding of ultracold collision processes between neutral atoms is crucial to the design and operation of 
atom traps, and to the development of novel quantum processes using trapped atoms Q . The elastic collision rate must 
be high enough to produce efficient thermalization during the evaporative cooling phase of magnetostatic trapping, 
whereas the inelastic collision rate must be small, because such collisions can generate energetic atoms and change the 
atomic state, hence destroying the trapping conditions and producing trap loss. Elastic collisions are also important in 
studies of Bose- Einstein condensates where they determine the mean field of the condensate [2( . Ultracold collisions 
are usually studied under weak trapping conditions in which the confining inhomogcncous magnetic field is cither 
ignored or is assumed to have a parabolic or harmonic spatial variation of sufficiently low frequency (typically 10 2 Hz) 
that it can be treated as uniform during the collision. However, recent interest in phenomena such as quantum phase 
transitions of 87 Rb atoms confined in three-dimensional optical lattices and far off-resonance three-dimensional 
optical lattices to store metastable argon atoms Q or to implement quantum logic gates and create highly entangled 
quantum states Q, involve conditions where the trapping frequency is typically 10 5 to 10 6 Hz and the tight trapping 
environment is expected to significantly modify the properties of the colliding system. 

In several existing calculations for tightly confined neutral atoms the exact interatomic potential is replaced by the 
regularised S- function pseudopotential 

V s (r) = ^aS(r)-^r, (1) 

where a is the scattering length and M is the reduced mass of the system. This potential reproduces the s-wave 
phase shifts in the Wigner threshold regime and also the correct asymptotic behavior of the wavefunction at large r. 
This enables an analytical solution for the case of a spherically symmetric harmonic trap to be obtained 6] , with the 
energy eigenvalues determined from the condition 

a - hf) - 1 1™ ( wE 4- n \ r( ^ + ^ } m 

l~ m = 2 tw \^ + Vn^fy (2) 



where uj is the trap frequency and £ = y^h/Moj is the effective range of the ground state wavefunction. The validity 
of this approach has been investigated by Tiesinga et al. y| for the Na and Cs systems by comparing the energy 
eigenvalues with those computed numerically using the best available full interatomic potentials. They find that 
this approximation is limited to sufficiently weak traps where £ S> a. For the case of two atoms interacting via a 
hard-sphere potential of range a, Block and Holthaus [8| have shown that the pseudopotential has the form with a 
given by J2J). Recently, two groups |9|, [lfj have advocated a model in which an energy-dependent effective scattering 
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length 



tan<$ (fc) 

acff(-E) = : (3) 



is introduced, where fik = y/2ME and 6o(k) is the s-wave phase shift for the untrapped atoms scattering in the 
interatomic potential. The energy eigenvalue condition 

^1 = m (4) 

is then solved self-consistently. The procedure given in J3J) and 10} reproduces the asymptotic wavefunction and the 
s-wave phase shifts even at energies above the Wigner threshold regime, and it is found that this self-consistent (SC) 
method works even when a G ff/£ 1. 

In this paper we analyze a system of two interacting atoms under strong harmonic confinement from the viewpoint 
of quantum defect theory (QDT) and derive a generalized SC model which validates |@J and extends it to non s-wave 
collisions. We then present two highly efficient computational methods for determining the bound state energies and 
eigenfunctions for the trapped atoms. The first method, based on QDT, is an adaptation of a technique developed 
by one of the authors (GP) for highly excited states in a modified Coulomb potential. The second method is a 
direct numerical solution of the radial Schrodingcr equation using a discrete variable representation (DVR) of the 
kinetic energy operator and a scaled radial coordinate grid. Both methods are applicable to a system of two arbitrary 
ultracold neutral atoms interacting in a harmonic trap and, as an example, we consider the case of spin-polarized 
metastable helium tightly confined in a spherically symmetric harmonic trap as this system has not been studied 
before under these conditions. 

Interest in collision processes in metastable 2 3 S helium (denoted by He*) has been generated by the quest to attain 
Bose-Einstein condensation in this system 0, ancl , subsequently, to understand these condensates 01 ■ Such 
condensates are novel in that they are the first excited-state condensates and open up new fields for investigation such 
as those of atomic correlations and the growth kinematics of the condensate. This experimental success has depended 
upon the correctness of the theoretical prediction [lj, Il5l Il6l Il7l 1 1 8| that the inelastic Penning ionization processes 
can be strongly suppressed through spin polarization of the He* system in the magnetostatic trap. 

This paper is organized as follows. In Sec. II the formalism is developed for collisions of two neutral atoms in an 
external three-dimensional isotropic harmonic trap, and the general nature of the energy eigenvalues and eigenfunctions 
of the resultant radial Schrodinger equation discussed. A quantum defect is introduced and shown to be analytic in 
energy. In Sec. Ill we formulate a generalization of the SC method, and the QDT and DVR computational methods 
are presented in Sec. IV. In Sec. V the QDT, DVR and SC methods are applied to ultracold metastable helium atoms 
under tight harmonic confinement. Results are obtained for s-wave collisions over a range of trapping frequencies and 
also for d-wave collisions in a 10 MHz trap. Finally, in Sec. VI, we summarize and discuss our results. 



II. TWO-ATOM COLLISIONS IN A HARMONIC TRAP 



Consider two atoms j = 1, 2 of mass Mj and position rj relative to the centre of the trap. For a central interatomic 
potential V(r), where r = \r\ = \r\ — r2|, and an isotropic harmonic trap, the two-atom Hamiltonian is separable into 
centre-of-mass and relative motions. The equation for the relative motion of angular momentum I is 



ip(r) — Eip(r) , 



(5) 



where E is the energy eigenvalue and the reduced mass M — M\Mij (M\ + M%). The trap potential VtrapM = 
Muj 2 r 2 /2 has been assumed to be independent of the atomic state, which is generally valid in far-detuned optical 
lattices As the interaction is spherically symmetric, ip( r ) has the form 



^(r) = -F kl (r)Y lr< 



(6) 



where Yi r 



are spherical harmonics and F k i (r) satisfies 



2M dr 2 



l(l + i)h 2 

2Mr 2 



-Mu) 2 r 2 



V(r) 



F kl (r) = EF kl (r) 



(7) 
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It is convenient to introduce the dimensionless variables p = r/£ and k = 2E/Tlui — £ 2 /c 2 and rewrite J7J as 

r d 2 1(1 + 1) 2 2V(p)' 



K — p 



F(p) = 0. (8) 



dp 2 p 2 huj 
For large values of p, V(p) is negligible and F(p) has the asymptotic form 

F a {p) = z l l 2+1 / 2 e-*/ 2 w(z)- z = p 2 (9) 
where w(z) is a linearly combination of the two independent solutions given by 

w\{z) — iF 1 (a; c; z) ; w 2 (z) — z C \F\ (1 + a — c; 2 — c; z) (10) 
on using the notation adopted by Luke [20|. In (|10|l 

i 3 K , 3 , . 

and i-Fi(a; c; z) is the confluent hypergeometric function 

^ r(a + n) r(c) z" 
iFi(a;c;2) = > — — — — — — -— . (12) 

^ T(a) T(c + n) n! 

For the case of the unperturbed oscillator the solution w\(z), which is regular at the origin, is bounded as z — > oo 
provided that a = —n r , where n r is a non- negative integer. The energy eigenvalues are then given by 

E = Ku(2n r + l+~); n r = 0,1,2,... , (13) 

where n r denotes the number of nodes in the corresponding radial wavefunction F^ (r) . 

In the presence of collisions, the energy eigenvalues E are no longer equal to Eq but can be written in the form 

E = hio{2n* r + l + ^). (14) 

Thus the effect of the collisions is to replace n r by 

n* =n r — (X (15) 

where p will be called a quantum defect and is in general not an integer; this quantum defect has previously been 
introduced by Blume and Greene (To|. 

The present problem, characterized by a long asymptotic trap region beyond the effective range of the interatomic 
potential, shows similar features to the modified Coulomb problem in which a long-range attractive Coulomb potential 
is supplemented by a short-range interaction. This is an extremely well-known situation in connection with the analysis 
of atomic energy levels and spectra where it is common to specify the energy of a level with quantum numbers nl in 
terms of an effective principal quantum number n*. The quantum defect is then defined to be the difference between 
n and n* . In this context there has been extensive development of a quantum defect theory, see Seaton [21| , and in the 
following, we demonstrate that equation (JSJ for the perturbed harmonic oscillator can be mathematically transformed 
into an equation that is identical in form to that for the modified Coulomb problem. This enables quantum defect 
theory to be directly applied here and valuable theoretical insight into the underlying behaviour of the energy level 
structure obtained. 

Asymptotically, the exponentially decaying eigenfunction @ is given by 

rfi c) -^(^ f) 

w 3 (z) = ip(a; c; z) = — — Wl (z) + w 2 (z) (16) 

1 (1 + a — c) 1 (a) 



where a = — n*, and as z — > oo, 



ip(a;c; z) — z a 2 F (a,l + a - c; --) , (17) 
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where 

„ . , T(a + n) T(c + n) z n 

2 F (a; c;z) = J2 -J^) T^T ^ ' (18) 

n=0 v ' w 

The quantum defect is an analytic function of energy for E > 0. To prove this we make the transformation 

F a {p) = z-V*Y(y) (19) 
where y — nz/8. The function Y(y) then satisfies 

d 2 A(A + 1) 2 



dy 2 y 2 y 



Y(y) = 0, (20) 



where A = 1/2 — 1/4 and 



1 K 

£ = -— ^; n* = -=n*+A + l= ^ — /i (21) 

with v = n r + A + 1. Equation (|20|) has exactly the form of the Coulomb equation for a bound state nl, in which n 
and Z have been replaced by v and A. The present problem differs from that studied by Seaton in that, here, A cannot 
be an integer. Consequently we introduce the two linearly independent solutions of 1|20|) 

(n*z) x+1 e~ z/2 (n*z) x+1 e~ z/2 
V) = ( r( ; A + 2) ^t(A + 1 - 2A + 2; z) = ( ^ + 2) m (22) 

and 

Y 2 {n*,X; y) = l " gj-^ 1 F 1 (-A - n*, -2A; z) = [U > w 2 , (23) 

where a = A + 1 — n* and c = 2A + 2 in (|10(1 . The functions Yi and Y2 are identical to the functions y\ and y 2 defined 
by Seaton and are analytic functions of e. We then introduce the general solution 



Y 3 (n*,\;y) = z x+1 e- z / 2 



a{n ) T{2xTY) Wl + ^ n) T^2Xj w ^ 



(24) 



where a(n*) and (3(n*) are analytic functions of n* (oc E), see (I14JI and l|21|l. On expressing w 2 in terms of wi and 
the exponentially decaying solution W3 in (|l(j[l . the condition for l|24|) to decay as z — > 00 is that the coefficient of w\ 
must vanish. This requires that 

o(n*) - A(n*,X) B(n*, A) (3{n*) = , (25) 



where 



for large values of n* , and 



1( - A, = r(n*%\t)^ - 1 + 0(1 " J| 



^A) = $f^| = (-!)' cot.M. (27) 
sin[7r(n* — AjJ 

In practice for the cases considered in this paper for which I < 2 , ^4(n*, A) remains very close to unity except for the 
lowest trap states. Therefore from <|25[) - (|27[l 

a(n*) sin(Tr^) - /3(n*)(-T) 1 cos(tt^) [1 + 0(l/n* 2 )] = , (28) 

and so fj, can be written in the form 

ri = a + bn* + cn* 2 + ... . (29) 

In general, the interatomic potential V(r) supports a number of bound states, rib say, for E < 0, so that the lowest 
trap state (E > 0) has n r = rib, that is, the lowest trap state eigenfunction has rib nodes. If however a pseudopotential 
is used such as Vs(r) in (Q, there are no bound states with E < and the wave function for the lowest state with 
E > has no nodes. In this case the number of nodes, n' r , and the quantum defect, fi' , are defined by 

n' r = n r — rib; ^ — H ~ n<b ■ (30) 

The differences between these two types of potential and the effects of their use are discussed in detail by Peach • 
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III. GENERALIZATION AND VALIDATION OF THE SC METHOD 

Beyond the effective range of the interatomic potential there is an extended region, typically 10 2 ao < r < 1 4 ao, 
where the effect of the trapping potential is extremely small and the asymptotic form of the wavefunction is essentially 
that of a free wave. Therefore in this region, to an excellent approximation ffcj(r) can be written as 

F k i(r) = kr[cosSi ji(kr) - sin Si raj (At)], (31) 

where ji(kr) and nj(fcr) are spherical Bessel functions, see [2^, and 5i{k) is the Z-wave phase shift. However in the 
previous section it has been shown that the correct asymptotic form of the wave function is obtained from (|19|l and 
1)240. so that (|31[) must correspond to 

F H (r) = C F a (p) = C z- 1 ^ [a(n*)Y x + (3(n*)Y 2 ] , (32) 

where C is a constant. In equation H32|) . the coefficients a(n*) and (3{n*) satisfy equations l|25|l - H27|) and so this form 
of the solution does decay exponentially at very large separations as it must for a bound state. We can ensure that 
the two forms ()31J) and (|32(l are identical throughout the region by examining their behavior for small values of r. 
The Bessel functions ji(kr) and ni(kr) in (|31|) take the form 

fcrjj(fcr) ^ (fcr) ,+1 ; krn^kr) ~ - ^ (fcr)~' , (33) 

whereas from IjlUf) and l|12fl . for small values of z = r 2 /£ 2 , 

1_ r(2A + 2) ' 2 " r(-2A) ■ 1 ' 



Since \ — 1/2 — 1/4, the behavior of i*fcj(r) as r — > in is 



' T" , i , ' 



e i+1 ra + |) r(-z + i) 

and then from H31|) . ll-i-il) and (|35|l we obtain 

1 /3(ra*) _ tan^j T(-l + \) Y(2l + l)r(2/ + 2) 



(35) 



(36) 



(n*) 2A +! a(n*) (^) 2/+1 r(f + §) [2T(Z + l)] 2 

where £ 2 fc 2 = 4n*. Combining this with the expression for /3(n*)/a(n*) derived from equations l|25|l - 127() . our final 
result is 



tan5;(fc) 
(C/c) 2 ' +1 

where n* = E/2huj and 



(37) 



For I = 0, (|38|l reduces to the result given by equations J3J and (0J. Thus not only have we obtained a generalization 
of the SC method to the case of / ^ 0, but we have proved mathematically that relation @ which was originally 
introduced empirically, is in fact rigorously true. 

The behavior of l|37|l for I > warrants some discussion. For potentials V(r) with the asymptotic form r~" at large 
r, the threshold behavior at small k of the phaseshifts 6i(k) is [23 

k 2l+1 ; n > 21 + 3 
tan<5j(fc)oc{ k 2l+1 \nk; n = 21 + 3 . (39) 
k n - 2 ■ n<2l + 3 



In our case, n = 6 and the behavior as k — > is given by 

tan 6i(k) 

tan<5;(fc) 



and 



fc 4 



-a;; 2 = 0,1 (40) 
W, l>2, (41) 



where a; and 6/ are constants. 
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IV. COMPUTATIONAL METHODS FOR TRAPPED ATOMS 



The accurate numerical determination of bound state energies and eigenfunctions for interacting atoms in a trap 
requires computational techniques optimized to deal with both the small inner interatomic region and the large 
outer trap region. Two such methods are now presented in some detail as they will form the basis of subsequent 
investigations into more complex processes involving trapped atoms. 

The first approach is an adaptation of a method (QDT) developed by one of us (GP) for highly excited states 
in a modified Coulomb potential. We consider the region r > ro(~ 20ao to 40ao) where the asymptotic solution 
specified in l(24*)l is only weakly perturbed by the presence of the interatomic potential V(p). We first of all generate 
this asymptotic solution throughout the region r > ro and then calculate a multiplicative correction factor in order 
to obtain an accurate solution of Equation (|2*0")) is written as 



dy 2 



g{y) 



Y(y) = 0; g{y) 



A(A + 1) 



(42) 



and the region is divided up into N ranges, y — y rll n = 0, 1, 2, . . .N, say. j/jv is chosen so that for y > yN the decaying 
solution Ys(y) can be evaluated by using its asymptotic form, see (|1 Tft and l)18[l. Within each range y n -i < V < Vn 
the following expansions are made 



oo / _ \ m Mo 

t — t m \ t — t 



(43) 



m=0 



and substituted into l|42l) . The coefficients a m can then be obtained by equating the coefficients of powers of (y — y m ) 
to zero, and the values of F(y n _i) and F'(y„_i) obtained from the range y n ~i <y<y n provide the input values for 
the solution in the next range y n -i < V < y n -i, etc. The choice of M in each interval clearly depends on its length, 
but it is found that the accuracy of the solution Y(y) is very insensitive to the precise choice of N and Mq. One of 
the y n is chosen to be at the outer turning point given by g(y) — 0, y = y a where 



Va = n 



A(A + 1) 



and yo and p a are defined by, sec lf^U|l and i|21fl. 

ya = n*(p f/2 ; y a = n*(p a ) 2 /2 . 



(44) 



(45) 



where p = r /^ and r a = t;p a . For y > y a (|42|) is solved to obtain Y 3 (y) only, but at y = y a the solution Y^{y) is 
introduced. We consider the Airy functions Ai(gr) and Bi(g) [23|, where y ~ y a and 



ig'(ya)} 1/3 (y-ya). 



(46) 



The functions Y^{y a ) and Y^(y a ) are approximately proportional to Ai(0) and Ai'(0), and a second solution Y±{y) of 
(|20|l . which is exponentially increasing as y — * oo, is chosen to be proportional to Bi(0) at y = y a . Hence 



and then the complex function 



- V3Y 3 (y a ) ; Yl(y a ) = -V3Y^y a ) . 



Y 5 (y)=Y 3 (y)-iY i (y) 



(47) 



(48) 



is propagated inwards over the range yo < y < y a . 

We now consider the evaluation of the correction factor R(x) where x = 1/p. The solution F a (p) of JSJ is written 

as 



F(p) = R(x)F a (p) , 



(49) 



where F(p) is given by (|19|l in which 



Y(y) = Y 5 (y), Po <p<p a ; Y(y) = Y 3 (y), p > p a . (50) 
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We also define the functions <p(p) and 4> a (p) by 



4>a(p) 



1 dF a (p) 



Fa(p) dp 



The function R(x) then satisfies the equation 



l d 2 R 
Ih! 2 



2x' 



dR 

dx 



2(j) a {r) x 



, dR 2V{x) 



dx 



R(x) = , 



which is solved using the boundary conditions 



0. 



(51) 



(52) 



(53) 



By making the choice for F a (p) given by l|50|) we ensure that 0o(/°) is a slowly varying function for all p > po and 
therefore that R(x) is also a slowly varying function of x over the range < x < xq, where xq = 1/po- Since r a may 
be of the order 10 3 ao to 10 4 ao, the range r < r < r a can be quite large. The choice of the complex function here 
is crucial because it behaves like exp(iz9) where $ varies approximately linearly with r thus guaranteeing that 4>a{p) 
changes slowly. Therefore the differential equation l|52|) for R{x) can always be solved very accurately throughout 
the outer region < x < xq , using a grid method containing a maximum of 66 points and finally the solution F^i (r) 
required is obtained by taking the real part of the function in Q49J1 . In the inner region, < r < ro, equation is 
integrated numerically outwards using the Numerov algorithm. By using an iterative process to pinpoint the precise 
value of the energy, the values of <fi(p) in l|51|l obtained from the two regions are matched at p = po. The value of ro is 
varied within the range 20ao < ro < 40ao and the results are shown to be extremely insensitive to its precise choice. 

The second approach to the solution of the energy eigenvalue equation (jSJ) is to use a scaled discrete variable 
representation (DVR) 26] (usually a Fourier grid) of the kinetic energy operator, and convert the bound state 
problem into one involving the diagonalization of a N x N matrix where N, for this single-channel problem, is equal 
to the number of grid points. Under a general real invertible transformation of the radial coordinate p given by 



t = u(p); p = u~ 1 (t) = U(t) 



equation JHJ becomes 



dt 2 



F(t) = nF(t) 



where 



and 



f(t) = 



dU 



' (j " /(*) 



Pit) 



Hi + i) +p2 + 2Vip 1 + nt) d 2 i 

p 2 Tko dt 2 



(54) 



(55) 



(56) 



(57) 



Equation l|55|l is solved with the boundary conditions F (ti) = and F(tiy) = by constructing a DVR using a set 
of basis functions {</> m (t)} and a finite set of coordinate points {t m } over the interval [£i, t^}. Then l|55|) is converted 
into N linear equations specified by 



53 [/ 2 (*i)^/ 2 fe) + P(tj)Sij] F(tj) = KF(ti) ; i = 1, 2, . . .N . 



For the Fourier basis defined by 



mir(t — ti) 



t'N — tl 



m= 1,2, ...N 



the DVR of the kinetic energy operator is 

J {-l)^[^c 2 {^{ l - J ))- C s C 2 {^{ l + ] ))} ■ i?j 
2(t N - tl ) 2 I i (2iv2 + 1) _ csc 2 (f) . i = j 



T 



(58) 



(59) 



(60) 
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The scaling transformation used 

t = u(p)=^y ; p^U(t)^Ct P (61) 

involves the parameters £ and p. Equation l|61|) gives ti and tjy in terms of and p^ where theoretically the 
boundaries are at p\ = and pat = oo. In practice, because of the extremely high potential barrier near the origin 
and the nature of the scaling, p\ is chosen so that p\£, = 2ao- The choice of pn depends on the number of states 
being investigated as the associated eigenfunctions must have exponentially decayed sufficiently on the outer boundary. 
Taking p^ = 15 allows us to calculate the positions of about the first 50 trap states accurately. The scaling parameters 
£ and p are chosen to ensure that a good proportion of the mesh points are inside the potential well around the point 
r = 7ao. On taking ( = 20ao and p — 10 about 17% of the scaled mesh points lie between r\ and C, and with N = 500, 
eight to ten digit convergence is obtained for the eigenvalues. 



V. APPLICATION: ULTRACOLD METASTABLE HELIUM ATOMS 



As an application of the computational methods adopted in this paper, we consider the case of spin-polarized 
metastable helium tightly confined in harmonic traps of various frequencies. The colliding atoms are in the 
molecular state for which we use the potential of Starck and Meyer |24|. This potential has a scattering length of 
156.777eto and supports 15 bound states. Calculated quantum defects for the 31 lowest trap states with I = are 
shown in Tables 1 and 2 for trapping frequencies ranging from 1 to 100 MHz. Also shown are the results obtained 
using the SC solution of (0J. Very recently Gadea et al. [25( have constructed a new 5 E S potential that supports the 
same number of states but has a scattering length of 291ao. Experiments ^lj, ^| suggest that these two scattering 
lengths represent lower and upper limits on the exact value. 

The QDT results are in excellent agreement with those obtained from the DVR method, with the absolute differences 
being O(10~ 6 ). Also the agreement between the results from the QDT and SC approaches is very good, the absolute 
differences being only O(10" 7 ) at 1 MHz, increasing to O(10" 5 ) for 10, 20 and 50 MHz, and O(10" 4 ) at 100 MHz. For 
the trapping frequencies considered, the bound states in the 5 E+ potential are relatively unaffected by the presence 
of the harmonic trap. The number of bound states with E < is still ni, — 15, but the most loosely bound state does 
show some unusual behavior, being shifted downwards at 1 MHz and upwards at 100 MHz. The trapping frequencies 
of 1, 5, 10, 50 and 100 MHz correspond to the ratios a/£ = 0.1167, 0.3691, 0.5221, 0.8255 and 1.167 respectively. 
The pseudopotential result based upon the use of an energy-independent scattering length a, breaks down at 
the higher trapping frequencies. Allowance for the variation of a e g(E) with energy E is essential in order to obtain 
the correct eigenvalues; a c g(E) changes sign as it passes through a divergence at k = 0.01356ciq , producing a rapid 
variation in a e g(E)/t; around E/hui « 165.6/^(MHz). 

As a test of the generalized self-consistent (CSC) method (|37|l we have considered the case of spin-polarized 
metastable helium atoms in a 10 MHz trap. We have chosen I = 2 since for identical atoms, parity considerations 
exclude odd values of I, and in this case the 5 S g potential supports 14 bound states and 62 = 2.383887 x 10 5 a 4 in 
i|4ip. Scaled energy eigenvalues n* calculated for the 31 lowest trap states using the QDT and GSC methods are 
given in Table 3. The two sets of results agree to six significant figures. The associated quantum defects p! for I = 2 
are much smaller than the I = quantum defects, increasing from O(10~ 4 ) for the lowest states to O(10~ 2 ) for the 
highest states considered. 



VI. SUMMARY AND DISCUSSION 



A QDT analysis of colliding ultracold atoms tightly confined in harmonic potentials has been undertaken and a gen- 
eralized SC model obtained. Two highly efficient computational methods are presented for calculating the bound state 
energies and eigenfunctions, one based upon QDT and the other on a scaled DVR method. The perturbed harmonic 
oscillator problem is characterized by a long asymptotic region beyond the effective range of the interatomic potential 
and the QDT method is very efficient for integrating inwards through this outer region. The radial Schrodinger 
equation for the relative motion of the harmonically confined atoms is transformed into that for a modified Coulomb 
potential and a quantum defect introduced that is an analytic function of energy. At large separations r, each eigen- 
function for the confined atoms is expressed as a product of the unperturbed harmonic oscillator eigenfunction and 
a residual function R(x) that is slowly varying in x = 1/r. The unperturbed function is calculated by dividing the 
region into a number of ranges and expanding the function in a power series within each range. The differential 
equation for R(x) is then solved accurately using a grid method requiring relatively few mesh points. The Schrodinger 
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equation is integrated directly through the inner region and the two solutions matched at an intermediate separation 
in the range 20 to 40 a . The eigenvalues arc determined very efficiently by iteration on an initial estimate obtained 
by extrapolating the quantum defects downwards in energy from highly excited trap states. 

The two computational methods have been applied to the case of trapped spin-polarized metastable helium atoms. 
Energy eigenvalues calculated with the QDT method agree closely with those computed directly from the radial 
Schrodingcr equation for the trapped atoms using a DVR method, and with those obtained using a self-consistent 
method (SC) involving an energy-dependent effective scattering length. 

The range of trapping frequencies and the number of trap states considered in the present investigation was moti- 
vated by the desire to test the validity and robustness of our computational methods and not by current experimental 
feasibility. The higher trap frequencies are beyond those presently used and the higher trap states may not be physi- 
cally realistic. The harmonic approximation to the confining potential of the optical lattice is only valid for the lower 
lying states of the atoms confined near the nodes or antinodes of the lattice. Also, the rate of quantum tunnelling 
from these higher states to neighboring wells in the lattice may be significant, thus requiring a calculation of the band 
structure type for the entire lattice. 

The present calculations need to be extended to include collisional loss processes and to study the loss rate as a 
function of trapping frequency. In particular, for trapped metastable helium atoms, it will be important to include 
the magnetic dipole-dipole interactions that couple the spin-polarized 5 £+ state to the 1 £+ state from which there is 
a high probability of loss through Penning and associative ionization at small interatomic separations. 
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TABLE I: Quantum defects [/ = fj, — 15 for the lowest states with E > and I — where n' r — n r — 15. Results are listed for 
1, 10 and 20 MHz harmonic traps calculated using quantum defect theory (QDT) and the self-consistent method (SC). 



— ; 

n 


1 MHz 




i n ft t r t 

10 MHZ 




20 MHz 






QDT 


SC 


QDT 


SC 


QDT 


SC 





— 0.06814498 


— 0.068145 


— 0.2202828 


— 0.220309 


— 0.3057509 


— 0.305860 


1 


—0.1002427 


-0.100243 


—0.3009027 


— 0.300927 


—0.3996139 


—0.399707 


2 


— 0.1237965 


— 0.123797 


—0.3540177 


— 0.354041 


—0.4593226 


—0.459905 


3 


— 0.1430405 


— 0.143040 


—0.3941769 


— 0.394199 


—0.5038833 


—0.503958 


4 


—0.1595646 


— 0.159564 


—0.4266773 


— 0.426698 


—0.5398083 


—0.539876 


5 


— 0.1741690 


— 0.174169 


—0.4540938 


— 0.454114 


—0.5701305 


—0.570193 


6 


— 0.1873241 


— 0.187324 


—0.4778824 


— 0.477901 


—0.5965138 


—0.596557 


7 


—0.1993344 


— 0.199334 


—0.4989491 


— 0.498965 


—0.6199703 


— 0.620024 


8 


—0.2104109 


— 0.210411 


—0.5178967 


— 0.517914 


—0.6411619 


— 0.641230 


9 


—0.2207071 


— 0.220707 


—0.5351471 


— 0.535164 


—0.6605457 


—0.660645 


10 


—0.2303387 


— 0.230339 


—0.5510070 


— 0.551023 


—0.6784511 


—0.678534 


11 


—0.2393959 


— 0.239396 


—0.5657064 


— 0.565722 


—0.6951230 


—0.695193 


12 


—0.2479504 


-0.247951 


—0.5794221 


-0.579437 


—0.7107489 


—0.710810 


13 


—0.2560604 


— 0.256060 


—0.5922930 


— 0.592307 


—0.7254758 


— 0.725529 


14 


—0.2637738 


— 0.263774 


—0.6044302 


— 0.604444 


—0.7394208 


—0.739466 


15 


—0.2711310 


— 0.271132 


—0.6159241 


-0.615936 


—0.7526786 


—0.752716 


16 


—0.2781659 


— 0.278166 


—0.6268491 


— 0.626862 


—0.7653274 


— 0.765366 


17 


—0.2849077 


— 0.284908 


—0.6372671 


— 0.637285 


—0.7774321 


—0.777473 


18 


—0.2913816 


— 0.291382 


—0.6472304 


— 0.647264 


—0.7890476 


—0.789086 


19 


—0.2976095 


— 0.297609 


—0.6567835 


— 0.656834 


—0.8002204 


—0.800253 


20 


—0.3036106 


— 0.303610 


—0.6659643 


— 0.666010 


—0.8109905 


—0.811021 


21 


-0.3094018 


-0.309402 


-0.6748057 


-0.674847 


-0.8213925 


-0.821424 


22 


-0.3149982 


-0.314999 


-0.6833365 


-0.683371 


-0.8314567 


-0.831488 


23 


-0.3204133 


-0.320413 


-0.6915817 


-0.691614 


-0.8412095 


-0.841236 


24 


-0.3256590 


-0.325659 


-0.6995634 


-0.699594 


-0.8506741 


-0.850702 


25 


-0.3307463 


-0.330746 


-0.7073012 


-0.707325 


-0.8598713 


-0.859898 


26 


-0.3356848 


-0.335684 


-0.7148123 


-0.714840 


-0.8688193 


-0.868843 


27 


-0.3404836 


-0.340483 


-0.7221124 


-0.722128 


-0.8775347 


-0.877560 


28 


-0.3451508 


-0.345150 


-0.7292154 


-0.729240 


-0.8860324 


-0.886055 


29 


-0.3496937 


-0.349693 


-0.7361338 


-0.736144 


-0.8943255 


-0.894348 


30 


-0.3541191 


-0.354118 


-0.7428791 


-0.742900 


-0.9024264 


-0.902449 
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TABLE II: Quantum defects // = fx — 15 for the lowest states with E > and I = where n' r = n r — 15. Results are listed for 
50 and 100 MHz harmonic traps calculated using quantum defect theory (QDT) and the self-consistent method (SC). 



n r 


QDT 


50 MHz 


SC 


QDT 


100 MHz 


SC 





-0 


.4487958 




-0 


.449380 


-0 


5764557 




-0 


.578201 


1 


-0 


.5537186 




-0 


.554150 


-0 


.6897921 




— 


.690958 


2 


-0 


.6208043 




— 


.621154 


-0 


.7650603 




-0 


.765956 


3 


-0 


.6716889 




-0 


.672001 


-0 


.8236159 




-0 


.824344 


4 


-0 


.7134089 




-0 


.713686 


-0 


.8724945 




-0 


.873108 


5 


-0 


.7491700 




-0 


.749404 


-0 


.9149564 




-0 


.915486 


6 


— 


.7807137 




-0 


.780921 


-0 


.9528010 




-0 


.953266 


7 


-0 


.8090982 




-0 


.809287 


-0 


.9871383 




-0 


.987552 


8 


-0 


.8350137 




— 


.835187 


-1 


.018702 




-1 


.01907 


9 


-0 


.8589405 




— 


.859100 


-1 


.048008 




— 1 


.04835 


10 


-0 


.8812266 




-0 


.881372 


— 1 


.075435 




-1 


.07574 


11 


— 


.9021318 




— 


.902268 


— 1 


.101266 




-1 


.10155 


12 


— 


.9218555 




-0 


.921982 


— 1 


.125723 




-1 


.12598 


13 


-0 


.9405553 




-0 


.940675 


-1 


.148980 




-1. 


.14922 


14 


-0 


.9583577 




-0 


.958469 


-1 


.171182 




-1 


.17141 


15 


-0 


.9753662 




— 


.975472 


-1 


.192445 




-1 


.19266 


16 


-0 


.9916661 




— 


.991766 


-1 


.212865 




— 1 


.21307 


17 


-1 


.007329 




-1 


.00742 


-1 


.232523 




-1 


.23271 


18 


-1 


.022416 




-1 


.02251 


-1 


.251490 




-1 


.25167 


19 


-1 


.036879 




-1 


.03706 


-1 


.269826 




-1 


.26700 


20 


-1 


.051062 




-1 


.05114 


-1 


.287583 




-1 


.28774 


21 


-1 


.064704 




-1 


.06478 


-1 


.304806 




-1 


.30496 


22 


-1 


.077940 




-1 


.07801 


-1 


.321536 




-1 


.32168 


23 


-1 


.090799 




-1 


.09087 


-1 


.337806 




-1 


.33794 


24 


-1 


.103309 




-1 


.10338 


-1 


.353649 




-1 


.35378 


25 


-1 


.115492 




-1 


.11556 


-1 


.369091 




-1 


.36922 


26 


-1 


.127370 




-1 


.12743 


-1 


.384160 




-1 


.38428 


27 


-1 


.138963 




-1 


.13902 


-1 


.398876 




-1 


.39899 


28 


-1 


.150287 




-1 


.15035 


-1 


.413262 




-1 


.41338 


29 


-1 


.161358 




-1 


.16142 


-1 


.427336 




-1 


.42744 


30 


-1 


.172191 




-1 


.17225 


-1 


.441113 




-1 


.44122 



TABLE III: Scaled energy eigenvalues n* for the lowest states with E > and I — 2 where n' r — n r — 14. Results are listed for 
a 10 MHz harmonic trap calculated using quantum defect theory (QDT) and the generalised self-consistent method (GSC). 



n r 


n' 

QDT 


GSC 


n r 


n 

QDT 


GSC 





1.749908 


1.749912 


16 


17.74052 


17.74052 


1 


2.749755 


2.749759 


17 


18.73948 


18.73949 


2 


3.749538 


3.749542 


18 


19.73840 


19.73841 


3 


4.749259 


4.749263 


19 


20.73728 


20.73728 


4 


5.748920 


5.748924 


20 


21.73610 


21.73610 


5 


6.748522 


6.748526 


21 


22.73488 


22.73489 


6 


7.748065 


7.748069 


22 


23.73362 


23.73363 


7 


8.747551 


8.747555 


23 


24.73232 


24.73232 


8 


9.746982 


9.746986 


24 


25.73097 


25.73097 


9 


10.74636 


10.74636 


25 


26.72958 


26.72959 


10 


11.74568 


11.74568 


26 


27.72816 


27.72816 


11 


12.74495 


12.74495 


27 


28.72669 


28.72669 


12 


13.74416 


13.74417 


28 


29.72518 


29.72519 


13 


14.74333 


14.74330 


29 


30.72364 


30.72365 


14 


15.74244 


15.74244 


30 


31.72207 


31.72207 


15 


16.74150 


16.74151 
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